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Abstract. 

We present new results from two codes, using tinite differencing and pseudo-spectral 
methods for the wave equations in (3+1) dimensions. We use a hyperboloidal transformation 
which allows direct access to null infinity and simplifies the control over characteristic speeds 
on Kerr-Schild backgrounds. We show that this method is ideal for attaching hyperboloidal 
slices or for adapting the numerical resolution in certain spacetime regions. As an example 
application, we study late-time Kerr tails of sub-dominant modes and obtain new insight into 
the splitting of decay rates. The involved conformal wave equation is freed of formally singular 
terms whose numerical evaluation might be problematically close to future null infinity. 



1. Introduction 

In the recent past, since the work of [1, 2, 3], hyperboloidal slicings have been used frequently 
for numerical computations in general relativity as an elegant solution to the outer boundary 
and the radiation extraction problem by including null infinity in the slicing [4]. The 
applications range from time-domain simulations of the scalar wave equation on Minkowski, 
Schwarzschild [3, 5] and Kerr background [6, 7] to compute quasi-normal mode oscillations, 
to study the late time behaviour of the solution, to study the influence of non-linear source 
terms, to compute the scalar self-force of point particles orbiting a Schwarzschild black hole 
[8] and to compute gravitational waveforms from extreme-mass-ratio inspirals [9, 10]. A 
new and interesting idea is to apply the framework to perturbation equations in the frequency 
domain [11]. Recently, the use of hyperboloidal slices even allowed to evolve the Einstein 
equations in axisymmetry including null infinity, see [12] and references therein. 

The spatial domain in the numerical solution of hyperbolic PDEs is typically truncated 
at a finite coordinate distance where artificial outer-boundary conditions are imposed and 
the outgoing radiation of the test field is extracted. This practice causes certain well-known 
conceptual and practical difficulties such as artificial reflections on the outer-boundary which 
can destroy relevant features of the solution. In this context, the use of hyperboloidal slicings 
offers an alternative by including the physical boundary on the numerical grid. 

Many numerical applications require a specific coordinate system in a compact domain 
and it is necessary to attach hyperboloidal coordinates at some transition point. In the original 
work of [1] this was accomplished by introducing a transition zone, which required the 
fine-tuning of many parameters for the chosen transition function and a disadvantageous 
hyperboloidal transformation which had a suitable asymptotic behaviour but caused the 
outgoing characteristic speed to drop in the transition region. This led to numerical problems 
in subsequent works |5, 6, 8]. In the recent work [10], appearing during the completion of 
this paper, the authors present a solution to the problem. 
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For this paper we developed a finite differencing (FD) and a pseudo-spectral (PS) 
code in (3+1) dimensions to test a new hyperboloidal transformation that simplifies the 
control over characteristic speeds for the Ist-order reduced wave equation on Kerr-Schild 
backgrounds. Thereby we can avoid the above mentioned problem by requiring the outgoing 
characteristic speed to be invariant under the compactifying hyperboloidal transformation. We 
demonstrate that by performing a numerical comparison with the attached hyperboloidal slices 
as originally used in [1]. As an example application we investigate the late-time decay rates in 
Kerr at the horizon, finite radii and null infinity and study the m -dependence of the splitting 
of certain sub-dominant modes which was found in the numerical studies of [6, 7]. More 
technically, we remove formally singular terms that appear in the conformal wave equation 
that might numerically be problematic to evaluate at the outer boundary J^+, where the 
conformal factor vanishes. Our framework is general enough to cover metrics of the Kerr- 
Schild form. Thus, our approach could be applied in the numerical analysis of quasi-normal 
modes in Vaidya [13], and more generally, in time-dependent non-axisymmetric Kerr-Schild 
spacetimes. 

The infra- structure and implementations we present here are comprehensive, allowing 
standard coordinates in the spacetime interior in (3+1), smooth matching of hyperboloidal 
slices and covering the large class of Kerr-Schild metrics with FD and PS techniques. A 
comparable work is [5] in which Minkowski and Schwarzschild backgroimds are treated 
on hyperboloidal slices in (3+1) with PS techniques using the aforementioned non-optimal 
transition zone. 

The paper is organised in the following way. In Sec. (2) we introduce our notation and 
convention of Kerr-Schild metrics and the scalar wave equation. At next Sec. (3) we briefly 
explain the hyperboloidal method as well as the particular hyperboloidal transformation we 
employed. Then we compare the resulting characteristic coordinate speeds and the scalar 
curvature for matched hyperboloidal slices with other methods. In section (4) we apply our 
method to solve the scalar wave equation with finite differencing, explain implementational 
details and compare different matching methods numerically in Sec. (5). We then proceed to 
explain the pseudo-spectral code and use this on a single hyperboloidal domain (no matching) 
to study polynomial decay rates of sub-dominant modes on Kerr. 

Note that we use a tilde as in r or g'^" to denote compatified hyperboloidal coordinates 
and components of tensors. Conformally rescaled tensors are denoted by a superscript $1 as 
in "Z'^ = l^^ /^l or '^g^^" = gf^^/Q"^. Derivatives wrt the Kerr-Schild radius r are denoted by a 
prime like dh/dr = h' and wrt the rescaled radius r by a dot dh/df = h. 

2. Wave equation in Kerr-Schild coordinates 

2.1. Kerr-Schild metrics and Kerr metric 

Kerr-Schild metrics are of the general form 

g^" = 7]^''' -2V{t,x^)Fr, (1) 

where 77'"^ = diag(— 1, 1, 1, 1) is the flat metric, {t, x*) are Kerr-Schild coordinates and l'^ is 
a null vector wrt r?'"' (^ g'^'lpiL = 0) of the form 

F = {-l,r), = {l,r)-2VF (2) 

and A;'* is the second null vector normalized to kf^li^ = 2{-^ l^PSij = 1). P and V fare free 
functions that characterise the Kerr-Schild metric in question. 

t Weassumeasymptoticflataesssuchthaty = M/r-(l + Vi/r + 0(l/r^))andZ*ni = l + n-li/r + 0{l/r'^), 
= x^r. In Kerr for example V = M/r{l + 0(l/r'^)) andn-l = 1 - 0^/2(1 - (n^)2)/r2 - ©(l/r*). 
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V{x') = M^ i^-V^n, (3) 



In Kerr with mass M and angular momentum aM we have 

.r^ , M 
s ^ 2 

where s := \/ (r^ — a^)^ + (2az)2, = ■^\/ r"^ — o? -\- s is the Boyer-Lindquist radius 
and r = + j/^ + the Kerr-Schild radius, and 



rBL^ + ' Tbl^ + ' Tb 



(4) 



Then is the future directed ingoing principal null vector and fc** the outgoing one. 
With our sign convention 

n-/ = + — , (5) 

r 

where n-l := n^l^^, = n* = x''/r is the radial unit normal. Note that is affinely 
parametrised, i.e. I'^W ^l^ = 0. Sometimes it can be useful to work with Z** := \/2W. Then 

2.2. Wave equation in Ist-orderform 

The hyperboloidal method, explained in the next section, includes a conformal rescaling of 
the metric. Therefore, we consider the conformal wave equation instead of = §. 
We reduce the conformal wave equation 

(□ - ^7^)^ = fl^^s^v. - r''(ff)VM -\ni, = Q (6) 

to Ist-order form by introducing the new variables ijjy := d^ilj, where r^(9) := g^^'^T^^ 
are the contracted Christoffel symbols and TZ the scalar curvature of g^" . In general 
F'* = (V^Z'' + (j)) holds for metrics of the form (1). In Kerr it can be shown that 

r^(5) = — (7) 

s 

This leads to the following system of evolution equations for the variables {V'i/=o,i,2,3, V'}- 
dotp = ipo, dotpj = djipo and 

1 

6' 

which is symmetric hyperboUc. For finite-differencing implementations of (8) the following 
flux-conservative form is more convenient for long-term stability of the numerical evolution 

-5°°5oV'o = 5°'a,^o + -^d,F' - \nil, where (9) 

F' :=v^(ff°Vo+ff'^V'j). (10) 

where finite differences are taken of V'o and instead of differentiating the evolution 
variables i/'o and i/'i directly. This, in turn, is more efficient if spectral methods are used 
to compute spatial derivatives and eq. (8) is preferred in that case. 

§ The scalar wave equation n^/) = is per se not conformally invariant, i.e. a rescaled solution '^i/> ;= ip/il is in 
general not a solution of the wave eq. of the rescaled metric ^g^^" := g^"^ /Q?, see for example appendix D of [14]. 
Note that 7?. = holds for Kerr-Schild metrics where P is geodesic like in Kerr or Vaidya spacetimes. 
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3. Hyperboloidal method 

Having explained the wave equation on Kerr-Schild backgrounds we are now ready to describe 
the hyperboloidal method in the following section based on [1]. We explain the steps 
necessary to transform a given spacetime metric g^^ in coordinates {t,x'^) and its t-slicing 
(t = const hypersiu'faces) to a conformal metric "gf^" in radially compactified hyperboloidal 
coordinates {i, x') and its hyperboloidal f-sUcing. It is often desirable to leave the coordinates 
{t, x*) untouched in the spacetime interior, inner domain, and attach the hyperboloidal domain 
smoothly at some transition point Ttr. The hyperboloidal method we employ consist of three 
transformations, a radial compactification and a conformal rescaUng of the metric 

= f2(r) with = 0, O ^ at (11) 
= with $~f2atJ^+ (12) 

where Q := dfl/dr and we assume and $ non-negative. For simplicity we set 17 = $ in the 
following, see Appendix B for Q ^ The conformal transformation extends the spacetime 
to include ^+ as its boundary [4] such that the conformal metric is regular there. For constant 
time slices to intersect we require the additional transformation 

I-2n-lV ,, I + 2n-lV 
t = t- h(r) with ^ <h' < K— (13) 

where the height function h{r) possesses a carefully chosen asymptotic singular behaviour to 
allow t-slices to penetrate J^"*", similar to time transformations from Schwarzschild slicing to 
horizon penetrating slicings, see for figure (2) in [15] for illustration purposes. The restrictions 
on h{r) in (13) are necessary for the hyperboloidal slices to be spacelike, where h' := dh/dr 

and / := ^1 + 2V{1 - n-l'^). 



3.1. Controlling the coordinate speed of in/outgoing characteristics 

In numerical applications of the hyperboloidal method, like for the hyperbolic system (8), it 
is often necessary to hmit the ratio ^ to obtain a stable numerical evolution, where Ar is the 
spatial grid spacing and the time step. This is the case, if the method of lines (MoL) and 
explicit time integration is used. Then the CFL-condition must hold c(r) ^ < v^^^, where 
c(r) is the coordinate speed of characteristics of the solution, explained in the following, 
and z/cFL the so called Courant number that is independent of the solution. Therefore, an 
efficient numerical computation is limited by max c(r) and computing time may be wasted in 
regions where c(r) < max c(r). On the other hand it may be useful to obtain more temporal 
resolution in certain regions, for example to resolve a particular feature of the solution or to 
increase the accuracy there. In the following we bring the height function derivative h'{r) in 
a particular form which allows to directly specify the outgoing characteristic speed c(r) as a 
function of radius in the hyperboloidal transformation, i.e. h'{r) — )• h'{r, c(r)) independently 
of ri. We also show how to use this feature for attaching hyperboloidal slices to an inner 
domain ||. 

The in/outgoing characteristic speeds in the direction of a hyperbolic system of the 
form dtU = A^diU + BU, where we assume U is vector- valued and A^,B are matrices, are 

II As mentioned in the introduction in the recent work [10] this problem is solved in a slightly different way. 
The authors do not explicitly consider the characteristic speeds but guarantee a smooth transition of inner and 
hyperboloidal coordinates through (bl — ^bl = i—f, see eq.(7) of [10], where tsL, r^^ are the Boyer-Lindquist time 
and tortoise radius. 
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given by two eigenvalues of A^rii. For the systems (8), (9) they are 

c^{9) ■■ = ± /i-9'°), (14) 

^ ±I-2Vn-l rd^i +1-2V 
1 + 2V ^ 1 + 2V ' 

where g"" = g^'^rii. We considered the limit n-l ^ 1 because in most parts of asymptotically 
flat spacetimes, see footnote 2, the two unit normals are almost aligned n-l ~ 1, e.g. 
in extremal Kerr for r > 1.2. Instead of proceeding with the straightforward but long 
algebraic operations to simplify c^{"g), we proceed in a more intuitive way by considering 
the characteristic speeds of the 1 st-order equation 

l^df^iP = (15) 

which applies to the characteristics of the PDEs (8), (9). Characteristics propagate on null 
rays of the background along which the solution is constant. The radial characteristic speeds 
of eq. (15) along the null vectors {—l'^), fc'* eq. (2) are 

cW=fcVfc° = n-?^^ cW^c+(ff), (16) 

c(0 = Fm/f = -n-l cC) = c-{g). (17) 

Under the compatifying hyperboloidal transformations (11), (13) c'^^^ , c^'^ change as 

cW=n.Z^^ c«^c+(5), (18) 

1 + 2V 1-h' \=^n-l ^ ' 

^^'^ =-"-^TX77^ cW^c-(a (19) 

1 + n'n-l 

where L := ^ = ^ From eq. (19) we can see that the conditions (11) on fi, tl make rS^^ 

vanish at (no ingoing characteristics). Eq. (18) shows that to make c*^''' non-vanishing at 
h' has to be of the form 

h'='±^(l-n-LH) ^'--^^^ '-^^il-n^LHir)), (20) 
1 - 2F ^ ^ 1 - 2V{r) ^ w; > V ^ 

where V(r) is V(x'^) along some arbitrary direction or some approximation. For practical 
purposes it is enough to consider the leading order term of V{x^), i.e. V{r) = 2M/r, as we 
do in the following, is a free function in the range from zero to one to modify the coordinate 
speeds in the hyperboloidal domain. The spacelike condition eq. (13) is guaranteed by h' of 
the above form, since LH > by definition ^. A simple choice might he H — 1. However, 
H aUows direct access to the characteristic speed c^'^^ = (1 - 2y)/(l + 2V) 1/H eq. (18) 
and it is more convenient to require 

Hi,)='^^-'±^. (21) 
1 + 2V{r) 1 - 2V{r) 

This condition is the core of our hyperboloidal method. It ensures that the outgoing 
characteristic speed is unchanged under the compatifying hyperboloidal transformation, see 

1 One may object that for h' of (20) in the limit CI ^ -^^2ni'^v ^ ' expanding this inequality in H at it 
becomes clear that this can only be if LH < 0. 
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Figure 1. Left: Characteristic speeds of tlie systems (8), (9) for KeiT a = 0.5 on t- 
slices (black dotted) r £ [2; 20] and compactified liyperboloidal i-slices witli f\j,^ = 20 
attached at rjR = 10 for the hyperboloidal matching condition (21) (gray) and for (23) with 
C = 20 of Refs. [1, 5, 6, 8] (light gray) for transition functions /-tanh-tan (light g. solid) with 
5 = l,s = 2 [5, 8] and /-exp (light g. dashed) [6, 1]. 

Right: ^'R,{f) for the conformal factors (transition functions f{f)) of Refs. [1, 5, 6, 8] (gray) 
and our choice (black, k = 6 dashed). The same five curves are shown in two light gray groups 
for shifted transition points Ttr = 2; 15. 



Fig. (1) (left) (gray curve overlapping black dotted curve), independent of the choice for fl. 
We set f2 in the following to be 

n{f)^l-f{r) with /(^~)=(^) e{X),X = P-f,„ (22) 

where fj^. is the transition point and w the width of the hyperboloidal domain. 

In the original work [1] and also adapted by [5, 6, 8] the height function derivative was 
set to 

V./H.^/ri + i^ + e^^^) (23, 

\ r ) 

where (7 is a parameter and / ranging from zero to one a transition function and the conformal 
factor was set to ^(r) = 1 — S — r\j!+. Hext has the correct asymptotic singular 

behaviour, where c+| j<r+ = S"^ /C"^ but the characteristic speed in the hyperboloidal domain 
is left uncontrolled in the transition zone and depends on the details of / (as the spacelike 
condition (13)). In this case c''^\ see eq. (18), drops like fi^ before h! in the denominator sets 
in to result in a non zero c'^^^ asymptotically. This effect is apparent in Fig. (1) (left) (light 
gray, light gray dashed), where we plotted (14) for our choice (20), (21) and for (23) as in 
Refs. [1,5,6, 8]. 



4. Wave equation on hyperboloidal slicing 

In the following we derive the exphcit expressions for "g'^", "7?. and r''( that appear in 
the conformal wave eq.(6) under the compatifying conformal hyperboloidal transformations 
(11), (12), (13). We identify all formally singular terms which could constitute potential 
obstacles for numerical evaluations at the outer boundary 



4.1. Transformation of the Kerr-Schild metric g^'^ and TZ, F'^ 

The metric components g'^^, g'^^ and the determinant of g^'^ transform in a straightforward 
manner, divergent terms cancel trivially, the expressions are given in Appendix A. The 
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regularity of 

V»^(i-.n.<v,(^)„..-._,..,.f^a4, 

at is guaranteed by the form of h! in eq. (20), since i^2A//r ^+ =• 1 + ^i^G* ^^-^ \ in 
asymptotically flat space times +and we obtain 



= (1 - 2n.;V)^-^ ( -HL ] {h' K^) . (25) 



The scalar curvature "7?. of '^' g^'^ is given by 
6 



- (^(1 - 2n-/V)7V + f2L(3- 2y) - f2Lf ''-fn-l^ , (26) 
"7^= - ^ (^(l-2n-ZV)7V + f2L(3-4y)) for Kerr, (27) 



1 

where "7 := 7/^! depends on the Kerr-Schild metric in question. In Kerr 7 ~ 2Af/s, 
compare with (7). In certain cases it is useful to distinguish the conformal factor $ and the 
rescaling factor Vl, then TZ changes as in Appendix B. Besides h! , the conformal factor has to 
be picked carefully to avoid steep gradients in "g^^'^ and "7?^ (which appear on the rhs of eqs. 
(8),(9)) as shown in Fig. (1) (right). 

Next we derive the expressions for F^'C^'g), where we first consider the conformal 
transformation behaviour of F''((7) under eq. (12): 

pM(n^) = "F^(5) + "C^, (28) 

where "C' is a tensor, "F'^(g) = T^'{g)/VL^, "n^ = (0,n*)/rj and = V/n. 

The transformation of F'^ ( " g) to compactified hyperboloidal coordinates involves the 
Hessian of (11),(13). We obtain 

r-("5) = r-("5)-^^V^ (29) 

5^ "9'-' = ^{-/^"(l - 2n./V) + 2h'-V{l - n-f) -2h'-, (30) 
{{Q."r - - 2n-fV) + n'{5 - 2^)) - 4f fl'n-lV} 

rli-l^ rl 1 1 

= {-L^h'i, - 2n.V) + 2.'. - n.P) - 2.'-, (31) 

~ 2n-;V) + Oi(5 - 2T/)) - AFtlLn-lV}. 

A closer look at the time component r°( "5) reveals that the third term in eq. (31) is 
divergent. It cancels with —2VlL"n^ — ~2ilL{—h' /il) that appears through the coordinate 
transformation "C^ — > "C^, see eq. (28). We remove these terms from "i^*' :— 
g^ug^x ri^i^A ^jjj ii^fi ^jj^ denote them as "H^ and "Cr, where "D^ is the remainder 
of the cancellation. 

F^( "5) = "r^ig) + "C^ - '^H^ + ''Z?^, (32) 

{2/i'L^, 0,0,0}. (33) 

r 

+ = 1 + + ^1) + <^(^-'') and in Kerr = 1 - *-^(a\n.f) + 0(r-4). 
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5. Numerical applications 

We are now ready for a numerical application of the methods described in the last section. 
We solve the system (9) using a (3+1) finite differencing (FD) code on attached hyperboloidal 
slices using g as in eq. (21) and h' as in eq. (23) as most commonly used in the literature and 
compare the evolution of the numerical errors. Then we solve the system (8) using a pseudo- 
spectral evolution scheme in (3+1), which allows higher accuracy than the FD code, to obtain 
new insights into the splitting of the late time tails of the solution [6, 7] and confirm known 
results about polynomial decay rates at finite radii and We also inspect the late-time part 
of the solution at the horizon, where oscillations have been predicted [16, 17] but until now 
not regarded in numerical studies. 

With both implementations we evolve non- stationary compactly supported initial data 
(ID) of the form ^(0, x') = 0, tpj{0, x') = and 

^o(0,i') = y'<""(i')e-(^-^'«)'/(2<^') and of the form (34) 

V'o(0,i') = e-(^*-*o)(2^-4)'5y/(2<^'), (35) 

where a,ro,XQ are parameters and y'''™(x') are the spherical harmonics wrt x\ We set 
the Kerr parameter M = 1 and define the local power index (LPI) to be LPI(['0]''") := 
dLogI [V']'™ l/dLogt, where [V']'™ is the Zm-spherical harmonic component of the field wrt 
Kerr-Schild coordinates. 

5.7. Finite differencing code 

Our finite differencing code is embedded in the Cactus parallelization framework [18]. 
We made use of the extension Llama [19] of the Carpet driver [20], which handles 
the domain decomposition of grids over processors and provides the required interpolation 
operations for boundary communication. The Llama code, comparable to the multiblock 
code of [21], allows to use a spherical grid with the "inflated cube" coordinates [22] in 
the angular directions, which consists of six overlapping coordinate patches each with two 
angular coordinates of the general form p = arctan(a;Va;'^), i ^ k. The finite differences 
are computed of "Vo and '^F^ in the radial and angular coordinates, then transformed to a 
common Cartesian basis and input to the rhs of eq. (9). The discretization of the fluxes "F* 
instead of the evolution variables tpi is beneficial for a long-term stable evolution. If used on 
a single Cartesian grid, see for example [23] and references therein, the form of eq. (9) is a 
necessary condition for the FD scheme to conserve a certain discrete energy norm [24]. Since 
we compute non-Cartesian FDs and use multiple overlapping coordinate patches, we require 
additional Kreiss-Oliger type artificial dissipation |25| to guarantees long-term stability. We 
use the method-of-lines with 4th-order Runge-Kutta time integration, 6th-order FD operators 
and add 5th-order artificial dissipation * (except at the radial boundaries) to the rhs of the 
evolution eq. for ipQ. 

5.1.1. FD code: Comparison of matching methods As a non-trivial test case we evolved off- 
centered Gaussian ID with Xq = (2.4, 1.2, 1.1), a = 1 on the dented f-slicing, see eq. (23), 
corresponding to the solid light gray curve in Fig. (1) (left) and on the smoother f-slicing, 
see eq. (21), corresponding to gray curve in Fig. (1) (left). We set the radial grid spacing to 
Ar = 0.05, the time step to A? = 0.04 and the number of angular grid points per patch to 
N^x Np = Alx 41. 

* This is different from the 6-patch-code of [21], where summation-by-parts I'D and dissipation operations are used 
at grid and patch boundaries. 
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Figure 2. Left: Evolution of -(/"o along x-axis for off-centered ID on dented hyperboloidal 
slicing (bottom) and smooth slicing (top) for Kerr a = 0.1. Slicing parameters and 
characteristic speeds as in Fig.(l) (left) light gray (/-tanh-tan) and gray curve (eq.(21)). The o. 
characteristic speed c+ drops in the hyperboloidal domain which is apparent in the solution. 
Right: Evolution of the errors of t/i along x-axis for /q = ID on dented slicing (squares) 
and smooth slicing (dots) for three resolutions Ai^rn,h'r = 0.0625,0.05,0.04 with fixed 
Np = 41. Inset: Convergence rate on smooth slicing. 
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Figure 3. Evolution of the components ['I'o]'""'^''* at for off-centered ID on smooth 
hyperboloidal slicing. Setup as in Fig. (2) (left top). The field components decay at late times 
with a power law ~ t ^" as apparent in the right plot. Sufficiently high resolution is necessary 
to see the correct power law at low amplitudes (dashed curves). 



In Fig. (2) the outward propagation of the initial Gaussian pulse along the i-axis is 
shown. The pulse passes the transition point r™ = 10 on the smooth slicing without apparent 
interferences (top panel). On the dented slicing (lower panel), the pulse slows down in 
agreement with c"^ of Fig. (1) (light gray) and leaves the numerical domain at later time. 
The larger gradients in the field at a given time, due to the bump in c+, cause the numerical 
error to be bigger as shown in Fig. (2) (right). 

The decomposition of the field ipQ at into spherical harmonics is shown in Fig. 
(3), an exponential oscillatory decay, followed by a polynomial decay with the asymptotic 
behaviour ^ t^" in agreement with the rule n = 1 + 2 for tp, i.e. n = ^ + 3 for ipo, 
derived analytically in [26] and confirmed in the recent numerical studies [6, 7]. We go into 
more detail about asymptotic decay rates in the next section, where we use pseudo-spectral 
methods and quad machine precision. 
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5.2. Spectral code 

To resolve the late-time part of the solution higher numerical accuracy and machine precision 
is required than we are able to obtain with the FD implementation described in the last 
section, see Fig. (3). For this reason we created a code that uses spectral methods to compute 
spatial derivatives and the method of lines with 4th order Runge-Kutta time integration. In 
the spectral domain the evolution variables are expanded into real spherical harmonics F'™ 
(angular direction) and Chebyshev polynomials C" (radial direction), i.e. into the basis of 
polynomials p'™" := y'™ C", where we use Gauss-Legendre,Gauss-Lobatto collocation 
points, respectively. The hyperboloidal method in combination with horizon penetrating 
coordinates does not require any boundary treatment and unlike the FD implementation we 
use hyperboloidal slices that cover the whole numerical grid, i.e. beginning at the first, ending 
at the last collocation point, where we use as for the FD code the polynomial conformal factor 
(22) with /c = 6, which is exactly represented in the coefficients space. If we wanted to 
use a hyperboloidal domain that covers the grid only partially, we would require two spectral 
domains that are joined at the transition point r™. Since the piecewise function (22) defined 
on a single spectral domain contains the Heaviside function which is badly represented in the 
coefficients space. 

The spherical harmonics and their derivatives Y^"^ .djY'-™ are expressed in Cartesian 
coordinates, similar to SpEC, see e.g. [27]. They are evaluated through another basis of 
harmonics ^'"^ = (n-^A/"]'™')', where A/"]'™' are constant complex null vectors AfiAfjS^^ = 
labeled by I and m, spanning the 2^ + 1 harmonics in each Z-eigenspace. We chose the same 
A/"]'™' as in [28] such that the F'"* and the are related by a Fourier transform in each 
Z-eigenspace 



ylrn _ ^Im ^Im' -im'mai ^'»" _ C_]^\"i / + m)\{l m)\ 

„^_, ' 47r(2l + l) ' ^ ' 

where a; = 27r/(2/ + 1). The derivatives 9j"I>'™ take the simple form 

9,$'" = a,n^(n'=AAf lNf'^\ (37) 

where dj-n'' = {6^^ — n^n^)/r and thus are tangential to the unit sphere. Given the 
decomposition of a function tjj into the polynomial basis p'™" we obtain the tangential 
component of the derivative " through c'jF''" C" and the radial component 9,. through 
y'^Sj-C" (or alternatively, through a recurrence relation for the Chebyshev coefficients). The 
complete partial derivative of ip is then given by 

dj^ = Uj drip + ^^ {dj)'ip. (38) 

Since we compute the rhs of eq. (8) in the physical space, we introduce an aliasing error which 

would cause high modes > /max to blow up during the evolution. For our implementation it 
was enough to erase the luAx-components of the rhs of i/' and at every step of the time 
integration to obtain a stable long term evolution and to leave the rhs of V'j untouched. 



5.2.7. Spectral Code: Application to Late-time Decay Rates in Kerr A generic non- 
stationary compactly supported perturbation on Schwarzschild background decays at late 
times as i^" , n = 2Z + 3, Price [29]. On Kerr background the picture is more complicated. 
Since Kerr is not spherically symmetric, neighbouring /-modes are coupled. Until recently 
there was a controversy in the literature about the value of n for Z > 4. In the simple picture the 
lowest Z-mode (compatible with the azimuthal and equatorial symmetries of the ID) generated 
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Figure 4. Left: Field components [i/)]'^"'^'* at ^+ for Iq ID in Kerr a = 0.5. The decay 
rates are in agreement witli tlie rule n = I + 2. 

Right: LPI of V around the horizon between = (1.8 - 2.21, 0, 0) for [Iq = 1, m = 1) 
ID in Kerr a = 0.1 with Nr = 90, Ng = 6 grid points. The asymptotic decay rates approach 
n = 5 in agreement with n = Iq + I + Z. 



by mode mixing during the evolution is dominating the solution at late times. Nevertheless, 
analytical work by Hod [26], Barack and Ori [16] predict a surprising dependence of n 
("memory effect") on the initially excited mode Iq, i.e. n — n{l, Iq). By several numerical 
studies [27, 30, 31, 32] it was justified that both pictures are correct and that the decay rate 
depends on the details of the initial value formulation, i.e. the particular coordinates (Kerr- 
Schild, Boyer-Lindquist) in which the spherical harmonics are defined and the particular 
choice of initial data. Therefore, we distinguish between BL harmonics Ibl and KS harmonics 
^Ks in the following through subscripts, where necessary. Another effect that was predicted by 
Barack and Ori [16, 17] are oscillations on the late-time solution at the horizon (OAH). 

Recently, the use of hyperboloidal slices in Kerr allowed Zenginoglu and Tiglio [6] to 
extend these numerical investigation of polynomial decay rates at finite radii to future null 
infinity (compactly supported non-stationary ID Iq < 5). Very recently, Racz and Toth [7] 
presented a more detailed study for various kinds of ID, including different fall-off properties 
towards null infinity up to /q < 6 where they investigated the behaviour of the sub-dominant 
modes as well (/-modes with decay faster than the slowest decaying mode). They could 
confirm numerically the m-independent rule, n = Iq + I + 3 for I > Iq and n = Iq + I + 1 
for I < Iq proposed earlier by numerical studies [31, 32] which covers the formulas derived 
analytically in [26, 16]. For the late-time tails at (TAS) the results of [6, 7] were found to 
fit the rule n — 1 + 2 for I > lo and n — lo for l<la — 2 which has been obtained analytically 
in [26]. The authors of [7, 6] found a radial splitting of the decay rates for certain values of Iq, 
I and m, where the decay exponent varies for observers near the black hole, distant observers 
and at Moreover, the effect was found to depend on the harmonic index m [7]. We 

observe such a splitting in our simulations as well (SPL). 

As an application of the hyperboloidal method we investigated the effects OAH, TAS, 
SPL on the late-time tails in Kerr. At first we check the convergence of the code for the non 
trivial test case of off-centered ID, see Fig. (6). 

TAS: Fig. (4) (left) confirms that the late-time tail decay rates at are in agreement 
with the rule n = I + 2, where the evolution of the 1 — 0,2,4 modes of ^ for = initial 
data is shown. 

OAH: In Fig. (4) (right) we show the LPI of tjj in the neighbourhood of the horizon 
between f = 1.8 — 2.2 (see figure). The field decays in agreement with the rule n = lf) + l + 3 
but no oscillations are present. The reason is the co-rotating azimuthal coordinate — 
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Figure 5. LPIs of [1/1]'==^ (left) and [V']'=^ '"=^ (right) for (Zq = 3, m = 0, 2) ID at 
radii between r = 2 — 12 (dashed lines) and r = 12 — 20 (solid lines). Shown is the splitting, 
see text, of the LPIs at different radii for m, = compared to m = 2, where no splitting 
appears. For A'^^ = 90, Ng = 7 grid points and a = 0.05. 
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Figure 6. Left: Convergence of L2-norm of Aip for off-centered ID with Nr = 

60, 65, 70, 75, 80, 85, 90, 95, where At/" is computed wrt A^r = 110. a = 0.9, Ng = 

10, At = 0.1 are fixed. Right: Convergence of constraint field ipi — ditp at Ng = 
13, 12,11, 10, 9 with fixed a = 0.9, Nr = 65, At = 0.1. 



4>BL — tBL^+ jlused in [16, 17], where the Kerr-Schild (/iks = arctany/a; as cj), which was used 
in [33], are independent of 

SPL: As we can see in Fig. (5) (left) the LPI of [^]'=3,™=o for (^^ = 3, m = 0) ID varies 
between n = 5 for 2 < f < 12 (dashed), n = 9 for 12 < f < 20 (solid) and n = 5 for f = 20 
at We observed the same splitting for {Iq = 3,m — 1) ID (not shown). This effect is 
not present for (^o = 3, m = 2) ID as shown in Fig. (5) (right) and for {Iq 3, m = 3) ID 
(not shown). The reason is that Kerr-Schild {Iq = 3, to = 0, 1) ID as well as the examined I 
mode contain the = 1 harmonic which decays with n = 1 + 1 + 3. The ratio of the ^bl — 1 
harmonic in the I = 3 Kerr-Schild mode decreases in the limits a — > and r — > oo, where 
we observe the field decaying with n = 3 + 3 + 3. It vanishes at where the BL and KS 
inclination coordinates agree and there the field decays like 7i = 3 + 2. By setting to > 1 we 
can exclude /bl = 1 from the initial data and the splitting vanishes. 

Jj According to [16, 17] the field decays at late times like ~ t^" e""^+"(''^'), where d ~ t is the Eddington- 
Finkelstein null coordinate and S1+ = a/{2Mr+), r+ = M + (M^ - a^y^^. 
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6. Conclusion 

We applied a compactifying conformal hyperboloidal (cch) transformation to the wave 
operator of Kerr-Schild metrics and derived all expressions, removed formally singular terms 
of the conformal wave equation for which the numerical evaluation at the outer boundary 
might be problematic. The hyperboloidal slices we used in that process let the outgoing 
characteristic speed invariant under the cch transformation independent of the Kerr-Schild 
metric in question. This was possible by analysing the characteristic speeds of the 1st- 
order reduced wave equation, whereby we obtained access to the arbitrary part H of the 
height function derivative for which the constraints on h' are automatically satisfied (suitable 
asymptotic behaviour, spacelike condition). This part allows to directly set the outgoing 
characteristic speed to the desired function of radius in the hyperboloidal domain eq. (21) 
and thereby to get more control on the efficiency and accuracy of the hyperboloidal method in 
numerical calculations, which we demonstrated through comparison with existing approaches 
in the Uterature. As an application we have numerically verified known decay rates of scalar 
test fields and their sub-dominant modes in Kerr at the horizon, finite radii and at and 
obtained new insights into the m-dependence of the radial spUtting of the late-time decay 
rates of certain harmonic modes. 

The asymptotics of hyperboloidal slices simplify the boundary treatment in numerical 
simulations and provide a clean solution to the wave extraction problem. Moreover, the cch 
coordinate transformation beyond the boundary can act like an adaptive numerical grid on the 
background, i.e. H and O can be used to uniformly spread grid points or act as a magnifier 
to increase the spatial or temporal resolution, say in the neighbourhood of a localised source. 
In addition to providing that, the expressions presented here could lead to applications of the 
hyperboloidal approach to Kerr-Schild metrics beyond Kerr. The next step in improving the 
PS code would be to implement two spectral domains to be able to use piecewise defined 
conformal factors. It would be interesting to test the stability of the FD and PS code with 
non-linear source terms as in [5] or to use a localised moving effective source [8] to compute 
the scalar self-force on a particle orbiting a Kerr black hole and finally to extend the code to 
handle gravitational perturbations on Kerr-Schild backgrounds. 
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Appendix A. Metric components ^g"^^, "5°^ and det "g^^ 

The cch transformation (11),(12),(13) change the metric components g'^^, g'-' of eq. (1) to 

n-Oj ^ n~Oj _2^yiO^~p ^ n~ij ^ n^ij _ nfi ^ j-, 

^rf^ = -h'Ln\ ""6'^ = 6'' + filL{2 + ftlL) n'n\ (A.2) 

[0 =l°-h'n-l, = r + n'f(lLn-l, (A3) 
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where = V/Q., in Kerr = "rBL/"s with "s{f,z) = n^s{f,z) = 

y(r2 _ (Oa)2)2 + (2a05)2 and "r^L = ^\/f2 _ ^2^2 + The determinant is 

detg^^ = l det"g^^ = l/L'^. (A.4) 



Appendix B. *7^ for $ ^ 1^ 



In certain cases it can be useful to distinguish the conformal factor <I> = JIA and the radial 
rescaling factor Q, i.e. A 7^ 1, say if the metric should be left unchanged but the numerical 
grid has to be adapted in some region. In that case nxgn^ gg becomes 



*7^= - 

*7e= - 



6 

A3 
6 

a¥ 



{l-2n-fV)N + L'l>{3-2V)j-L^n-l"j^ , (B.l) 
{1 - 2n-l^V)N + L^{3 - 4V)\ in Kerr, (B.2) 



where TV" := ^{^"r - = L^{Li^ + 0$)f - L$, Li := 20 + QSlfL. The conformal 
transformation r^( "g) T'^i^'^g) of eq. (32) is given by 

PM("A^) = "f''(g)/A2 + "^4^ - ''H^/X^ + ''^Di", (B.3) 

"^C" = -2$L(''^n'^ - 2n-/ '"y /^)/A2, (B.4) 

^^i)" = {2/i'L^ (f A/A + 1) , 0, 0, Oj/A^. (B.5) 
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